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SUMMARY 


This report documents a portion of the work accomplished in Contract 
NAS5-11385 to investigate the use of optimal multi-impulse trajectories as a 
nominal about which one may expand using asymptotic expansion techniques to 
obtain approximations to optimal low thrust trajectories. The work documented 
herein consists of the analysis and description of an optimal 3-impulse trajectory 
program. Unlike other optimum multi-impulse programs, this one incorporates 
a patched-conic trajectory model and is specifically designed for compatibility 
with the subsequent addition of the low thrust expansion approximation. A report 
presenting the equations for the expansion of an impulse into an optimal finite 
burn is being published concurrently with this report. 
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I. introduction 

In this report we describe in detail a procedure for obtaining an optimal 
three-impulse, patched conic trajectory, from orbital motion about a 'launch" 
planet (taken as the earth) to orbital motion about a target planet. A computer 
program implementing this procedure has been constructed, and some numerical 
results are presented and discussed in the last section. 

Since we are using a patched conic formulation, the trajectory will con- 
sist of a sequence of Kepler arcs along each of which the vector A, adjoint to 
the velocity vector, and its time derivative A propagate as solutions of the 
Kepler variational equations. The end points of the successive Kepler arcs de- 
fine a sequence of five instants in time at which various boundary and transversa- 
lity conditions must be met: 

1. At t = t insertion is made from orbital motion about the earth 

o 

to an escape hyperbola. It is assumed that this insertion will be made by a high 
thrust, low specific impulse engine, with an option to jettison after use. 

2. t is the time of exit from the earth's SOI (sphere of influence). 
Between t Q and t the motion of the vehicle is along an earth focused Kepler 
hyperbola. At the SOI the attracting gravitational field switches from that of 

the earth to that of the sun, and at the same time transformation from geocentric 
to heliocentric coordinates is made. The vector A and all state variables are 
continuous across the SOI, but the transversality conditions require a discontinuity 
in A at this point. 

3. tj is the time at which a midcourse impulse is executed. Between 

t and t the vehicle moves along a sun focused ellipse. It is assumed that this 
X I 

midcourse impulse is made by a low thrust, high specific impulse engine, again 


- 1 - 



with an option to jettison after use. The presence of the impulse implies a dis- 
continuity in velocity at t ; the transversality conditions require continuity for 
• ^ • 

both A and A and also that A and A be perpendicular. 

4. t is the time of entry into the target’s SOI. Between t T and 

N I 

t^ the vehicle moves along a second sun focused Kepler ellipse. At t^ the 
attractive gravitational field switches from that of the sun to that of the target 
planet, and a transformation from heliocentric to planetocentric coordinates 
is made. All state variables and A are continuous, but, again, the transversality 
conditions require a discontinuity in A. 

5. t^ is the time for the retro-maneuver. The vehicle moves along 

a target-focused hyperbola between t and t., and is inserted into orbital 

N f 

motion about the target at t^. It is assumed that this maneuver is executed by 
a high thrust, low specific impulse engine, with an option to jettison. 

The three engines used for the impulses at t Q , t^ and t^ are characterized 
by their exhaust velocities c c and c respectively. The dead weights 

1 Z O 

(which may be jettisoned) are taken proportional to the initial mass, with propor- 
tionality factors k k , k respectively. Assuming that all three jettison 

X Z o 

options are exercised (and this is the case for the numerical results) the ratio of 
the final mass to the initial mass is easily shown to be 

II [l - (1 + k.)(l - e ^V C i) ] 
i=l 

where the Av., i = 1, 2, 3, are the magnitudes of the velocity discontinuities 
associated with the impulses at t Q , and t^ respectively. The goal of the 
optimization procedure will be to maximize this mass ratio (or minimize its 
negative), subject to a set of initial and terminal constraints on the earth and 
target orbital motions, by proper selection of the three impulses. 


( 1 . 1 ) 
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Because of the constraints and discontinuities associated with the five 

special points, defined by t , t„, t„ t and t_, and the fact that three co- 

O X I N I 

ordinate systems are involved a rather elaborate system of notation has been 
developed and is given in Section II. Also included in Section II is the pertinent 
information on the coordinate transformations used, the equations of motion for 
the vehicle, the differential equations for the adjoint variables and the transversa- 
lity conditions (without derivation) which are explicitly used in Section III. 
Derivations of the transversality conditions is made in Section VI, and a more 
complete discussion is given there. 

The procedure for constructing the optimal three impulse orbit to orbit 
trajectory consists of a double iteration in which the time tj and the position 
Rj of the midcourse impulse are the control parameters to be differentially 
corrected at the end of each pass through the double iteration. We outline this 
procedure, discuss the role of terminal constraints and motivate our choice of 
control parameters in Section III. An initial guess on the control parameters is, 
of course, required to start the double iteration. The method used for this 
initialization is presented in Section IV. 

The first part of the double iteration is an iteration (also requiring Initia- 
lization, given in Section IV) to find a sequence of Kepler arcs satisfying the ter- 
minal constraints and the current values of the control parameters. This iteration 
is described in Section V, where also complete details are given for the choice of 
terminal constraints implemented in our computer program. The second part of 
the double Iteration involves the determination of inconsistencies in the propaga- 
tion of the adjoint variables due to non-optimality of the current sequence of Kepler 
arcs. In our formulation these, inconsistencies appear as failure of the adjoint 
variables to satisfy the transversality conditions at t^. The calculations involved 
In the second part of the iteration are given in Section VI. The last part of the 
double iteration procedure is to differentially correct the control parameters, and 
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the basic idea involved is given in Section III. Since the optimization problem 
is sensitive, and the procedure for obtaining initial guesses is somewhat crude, 
a sophisticated differential correction procedure is required. The one selected 
for our program is described in Reference 1. 

Finally, in Section VII we present and discuss the numerical results for 
an optimized three impulse trajectory from an orbit about the earth to an orbit 
about Mars. 



II. 


NOTATION AND BASIC EQUATIONS 


We use a number of subscripts and superscripts, listed below, which 
may be appended either individually or in combination to the various symbols 
which follow: 

Subscripts X, I, N refer to quantities evaluated at exit from the 
earth's sphere of influence, the midcourse impulse point and entry to the tar- 
get's sphere of influence, respectively. 

Subscripts o, f refer to quantities evaluated at the initial and final 
times respectively. 

Subscripts E and T refer to quantities associated with earth and target 
respectively, including references to the earth centered and target centered co- 
ordinate systems. No subscript is used for solar reference. 

Superscripts + and - refer to evaluation just before and just after a 
discontinuity. 

R, r 

R, v 

t 

A, X 

A 


are the position vector of the vehicle and its magnitude, 
except that R^ also refers to the position of the midcourse 
impulse. 

are the velocity vector of the vehicle and its magnitude, 
is time 

* 

the primer vector, adjoint to R, and its magnitude. A 
subscript V is used to distinguish this variable in 
Section VI. 

time derivative of the primer vector and the negative of A , 

R 

the vector adjoint to R . 
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radii of the earth and target spheres of influence, res- 
pectively. 


m e . /*. 


G times mass of earth, sun and target, respectively. 


• •« 
P, P, P 



AV 1* AV 2 ’ 


planetary position, velocity and acceleration vectors. 

unit polar vectors for earth and target. 

exhaust velocities of the engines used to produce impulses 
at t^, tj and t^, respectively. 

masses of these engines and associated propellant tankage 

divided by m . 

J o 

AV magnitudes of velocity differences produced by these 

O 

impulses. 


n 

a, e, i 
H, h 


cost function = - final mass/initial mass, 
semimajor axis, eccentricity, inclination, 
angular momentum vector and its magnitude. 


a 


longitude of ascending node of planetocentric hyperbola on 
planet's equator. 


to angle from ascending node to spacecraft on planetocentric 

hyperbolic trajectory, measures in the direction of motion. 

y planetocentric flight path angle. 

Subscripts E and T refer to orbital motion about earth and target, res- 
pectively. An additional subscript H refers to the hyperbolic trajectories. 


Other symbols, used locally, are introduced as needed. 


The planetary coordinates and velocities will be needed in the coordinate 
transformations and are available in various ephemerides. In our program we 
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( 2 ) 

used an approximate analytic ephemeris from J. P. L. v . The coordinate sys- 
tems used in our program for the earth focused, sun focused and target focused 
segments of the trajectory are, respectively, geocentric, heliocentric and tar- 
get centered, all with their axes parallel to a geocentric equatorial system. To 
transform a vector V and its time derivative V among these coordinate systems 
at a time t , we use the formulas 

V <t ) = V E (t ) + P E (t ) = V T (t ) + P T (t ) (2. 1) 

V (t) = V E (t) + P E (t) = V T (f) 4- P T (t) 


The equations of motion of the vehicle are 


V-" E V r ! 

t ^t^t 
o X 

R r 3 

‘ X StSt N 

**T = _/i T R T /r T 

V tst f 

in accordance with Eqs. (2. 1), 

R x r ex + %x- 

A x ■ *EX + ^EX 

r n _ r tn + p tn’ 

^N _ ®TN + ^TN 


( 2 . 2 ) 


(2.3) 


implying continuity of the state variables, R, R over the spheres of influence. 
The accelerations are, however, not continuous. Us ing the approximate Keplertan 
nature of the motion of the earth and target about the sun, we have 
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EX 


(m + p e ) 


EX 


EX 


(2.4) 


TN 


(M + M t ) 


TN 


TN 


where the p's are the magnitudes of the P's, and hence 


EX 


AR„ - R„ - b" = -yi + - V s " + <^E> ^ 

r X r EX EX 


X X X 


R. 

• • • ♦ -f- • • ^ 

AR N = R N- R N =+(i — 

r N 


R 


- M, 


TN 


T 3 
r 

TN 


TN 

(M^ T ) — 

P TN 


(2.5) 


The primer vector A satisfies the Kepler variational equations: 

A - - m e A/r* + 3 4 e (A- R e ) R/r* t # *t^ t x 


A = - n A/r 3 + 3 [i (A*R) R/r 5 


t ^t 
X N 


( 2 . 6 ) 


A = - H T A/r 3 + 3 m t (A' R t ) R T /r^, ^tstj 


Transversality conditions at the sphere of influence require discontinuities in 
A at these points, given by 


W R EX <V Ah X )/(R EX- R EX» 


; + 


A N A N + R TN (A N ' AR N ) ^ ^TN ' R TN J 


(2.7) 
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The primer vector A is continuous at the spheres of influence. Transversality 
conditions at the midcourse impulse point require 


• 4 - • - 

r i“ r i 

A= A = A continuous 

I I Av„ 


A = A 


continuous 


( 2 . 8 ) 


WAjAj-0 


The initial and final constraints provide the boundary conditions for the 
state variables, while the initial and final transversality conditions provide 
boundary values for the primer vector and its time derivative. Discussion of 
the state variable constraints appears in Sections IV and V and the initial and 
final transversality conditions are given in Section VI. 
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in. THE DOUBLE ITERATION PROCEDURE 


The optimization of the patched conic trajectory of a space vehicle, from 
orbit about the earth to orbit about a target planet, is a two point boundary value 
problem involving the six state variables, R and R , and their six adjoint vari- 
ables, A and A as dependent variables, with time as the independent variable. 
There are, thus, thirteen parameters associated with each terminus of the tra- 
jectory: the terminal time and the terminal values of the twelve dependent vari- 
ables for a total of twenty six terminal values. In this problem t is a signi- 

o 

ficant parameter since the time origin determines the configuration of the sun 
and planets, in whose gravitational fields the vehicle moves. Since there are six 
state variables a maximum of seven constraints (six state and one time) may be 
specified at each of the terminal points. For each unspecified , or "free", func- 
tion of state variables and time at either terminus a transversality condition on 
the adjoint variables holds. Hence, in all cases, there will be seven initial and 
seven final conditions on the aggregate of state variables, adjoint variables and 
time. 

The first basic idea involved in the double iteration procedure is to de- 
couple the state and adjoint variables in such a way that we can find a sequence 
of Kepler arcs satisfying the terminal state constraints and then generate a 
time history of the adjoint variables along these Kepler arcs. This idea led us to 
characterize the midcourse impulse by the time t^ and the position R^ at 
which it occurs. These parameters t^ and R^ are regarded as "free", that 
is as unconstrained, although we must have current estimates for their values 
at each entry into this calculation. The second basic idea is now to decouple the 
trajectory problem (finding the sequence of Kepler arcs) into a "forward" and a 
"backward" part: 
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(1) Find Kepler arcs from initial state and time values to t^, R^. 

(2) Find Kepler arcs from t , R^ to final state and time values. 

The solution to each of these problems requires a total of four terminal 
values consistent with the terminal constraints imposed on the problem, and 
augmented, if necessary, by current estimates of "free" terminal parameters. 

In our implementation of this procedure the initial and final times are specified, 
together with three initial state constraints and three final state constraints. 

There is, therefore, no need to obtain current estimates on "free" terminal para- 
meters. The modifications necessary for other distributions of terminal con- 
straints are discussed briefly at the end of this section. 

In the next section we give the procedure for obtaining an initial estimate 
for tj and R^, as well as the initialization procedure required for the trajectory 
iteration described in Section V. An iteration procedure is used for the trajectory 
calculation because the changes in coordinate system and the equations of motion 
across the earth and target spheres of influence render the problem awkward 
for analytic treatment. 

» 

The details for the generation of time histories for A and A are given in 
Section VI. Since the sequence of Kepler arcs obtained from the trajectory calcu- 
lation cannot be assumed optimal, one would not expect to be able to propagate 
the adjoint variables by the Kepler variational equations along these arcs and also 
match all the transversality conditions. In fact, one would expect that exactly 
four transversality conditions would be violated because the eight time and ter- 
minal state constraints assumed would (in general) yield a unique sequence of 
Kepler arcs for each choice of the four free parameters tj, R^. Further, 
since these are the four free parameters for the trajectory calculation we choose 
the transversality conditions at t^ 
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Aj continuous 


j A • A 

> A .i- -Vr 5 - 


(3.1) 


to measure the "mismatch". After completing a pass through both the trajectory 


and A, A calculations, differential corrections for and R^ are computed 
to give the current estimates for the next pass. The principle involved is simple: 
defining the 4-vectors X and U as 

■'N 


r. + 

A i - A i 


x = 


*V(A+ +Aj)A I 




(3.2) 


r ^ 
4 


j 


u = 


R 


V J 

one must find U such that X vanishes. Since X is a very complicated non- 
linear function of U, an analytic solution is not feasible. Linearized corrections 
to U may be obtained from 

r ^ 




AU = 


AR. 

j 

V J 


a x 


a u 


V J 


-l 


x 


(3.3) 


and we update t and R^ by 


New tj = old t + Atj. 
New R t = old Rj + ARj 


(3.4) 
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The success of this method depends on having current estimates to t^ and R^ 
that are close to the values for which X vanishes. If this condition is not satis- 
fied a more sophisticated correction procedure is required. The details of the 
correction method used in our program are given in Reference 1. 

Sufficient detail has now been given to enable us to summarize the double 
iteration procedure succinctly, as follows: 

I. Generate, by iteration, a sequence of Kepler arcs satisfying 

initial and final time constraints and three constraints each on the initial and 

/ 

final states, and passing through the current estimate of R^ at the current 
estimate of t^. 

II. Generate a time history of the adjoint variables, A and A, along 
the Kepler arcs generated in I, in such a way as to satisfy the initial and final 
transversality conditions and the transversal ity conditions at the earth and target 
spheres of Influence. 

III. Compute the X vector, which measures the failure of the adjoint 
variables to satisfy the transversality conditions at t^, R^ and then differentially 
correct tj and R^ to drive X to zero. Return to step I and iterate until con- 
vergence is obtained. 

As noted above this procedure refers to the case for which the initial and 
final times and three state constraints at each of the terminal points are satisfied. 
To generalize the procedure is not difficult. For example, to leave one or more 
of these eight constraints unspecified requires that we impose the corresponding 
transversality conditions on the A, A calculation, and augment the "free” con- 
trol variables with the unspecified parameters. That is, we augment the X 
vector, using the additional transversality conditions and the U vector by the 
additional free variables. An initial estimate for the additional free variables 
must be made, with subsequent current estimates coming from Step III. Just 
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how the initial estimate is made would depend on the selection of the free variables. 
To impose more than four constraints at either terminus, one selects four for use 
in the trajectory calculation and uses the mismatch in the remaining ones to aug- 
ment the X vector. Note that if more than a total of eight terminal constraints are 
imposed the four free parameters t and will, in general, be insufficient to 
eliminate all mismatches. 


- 14 - 



IV. INITIALIZATION PROCEDURES 


Because apse and nodal crossings are known to represent optimum im- 
pulse points for certain types of maneuvers, procedures for obtaining initial 
estimates of tj and R^, given t Q and t^, for each of these points were in- 
corporated into the program. For impulses at either of the apses a heliocentric 
two-body trajectory between the ecliptic projections of the launch and target 
planets is computed by solving Lambert's problem. The indicated apse (i. e. , 
either periapse or apoapse) time and position are then employed as the guesses 
for tj and if that apsis is passed between the specified times t Q and t^. 

The nodal crossing points are obtained by forming the cross product of the 
angular momenta of the two planetary orbits. The intersection of the indicated 
node with the heliocentric two-body trajectory resulting from the solution of 
Lambert's problem above then yields the estimates of time and position t^ and 
Rj. If the indicated apse or node is not contained within the ecliptic two-body 
trajectory, a search for one of the other possible choices is then made and, if 
available within the trajectiry, is used to commence the iteration. 

Another initialization procedure must be carried out at each entry to 

Step I of the double iteration procedure. As will be seen in the next section, 

• * 

the trajectory iteration requires current estimates of R and R , the 

velocities of the vehicle at exit from the earth's SOI and at entry to the target's 

SOI, respectively. In obtaining our initial estimate for these velocities, we 

assume that the initial estimates for t^ and R^ are already available. Using 

again the ephemeris values for P and P (positions of earth and target at 

rt 

t and t, respectively), we solve two heliocentric Lambert problems : 
o I 

P„ to R t in transit time (t T - t ) 

Eo I ' I o' 

( 4 . 1 ) 

Rj to P,^ in transit time (t^ - tj) 
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to obtain initial and final velocity vectors P , P (which, of course, differ 
from the planetary velocities at these times, and hence are denoted by bars). 
These Lambert problems give a gross approximation to transfer from the earth 
to the target planet via at time tj, and have two characteristics that are 
of use to us : 


(1) the composite trajectory pierces both spheres of influence; mani- 
pulation of velocities at t , R^ could easily miss the earth and target altogether. 

(2) it turns out that the and associated with these Lambert 
problems are excellent approximations to their final values after convergence of 
the iteration. 


Thus our first estimates on t , R , t , R are obtained simply by 

A A JN JN « 

propagating one heliocentric Kepler ellipse forward from P , P to a time 

Eo Eo 

t such that 

.X. 


I R X- 


Eo 1 


= r 


E 


(4.2) 


and a second Kepler ellipse backward from P . P„ to a time t such that 

TF TF N 


l R N “ P TF 


• • 
and evaluating R^. as the velocity on the first ellipse at t and R^ as the 

velocity on the second ellipse at t . 

N 

The information so far generated in this section, t , R , t , R , t , 
9 I I X ' X N 

R n provides the input for the first pass through Step I of the double iteration 

procedure. Updated values of t ,, t^ T , and R. T , for subsequent passes, 

X N X N 

are calculated just before leaving Step I, and updated values of T^ and R^ are 
generated in Step III. 


(4.3) 
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One further Initialization is required, and should, perhaps, have been 
mentioned first. The coordinate systems used in the analysis depend on t and 
tj, and are defined in terms of the relations among the earth-equatorial, ecliptic 
and target-equatorial systems at these times. Since many terminal constraints 
are given most conveniently in terms of elliptic orbital parameters (e. g. inclina- 
tion) relative to the planets involved, all such coordinate systems must be ini- 
tialized. Note that if the generalization is made to "free" t and/or t., this 

o f 

initialization must precede every entry into Step I where the terminal constraints 
are used in the trajectory calculation. 
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V. 


THE TRAJECTORY ITERATION 


The trajectory calculation presented in this section is a modified version 
( 3 ) 

of that given by S. Pines' . This report deals with a direct transfer from earth 
orbit to Mars orbit with no midcourse impulse, and uses terminal constraints 
which are different from ours. The basic idea of using positions and velocities 
at the spheres of influence to control the iteration is retained. 

As noted at the end of the previous section, the input quantities for the 

trajectory iteration are t , t , t T , t , t,, R v , R. T and R T , of which t 

oXINfXN I o 

and t^ are assumed given and fixed. The remaining parameters are estimated 
by the initialization procedure for the first pass through the trajectory iteration. 

For succeeding passes t , t , R and R are updated at the end of the tra- 

IN A JN 

jectory calculation, while t^ and R^ are updated in Step III of the double itera- 
tion procedure. 

Before presenting the iteration, we specify the constraints to be imposed 
at the terminal times and t^. The constraints selected for our program are 
the inclination and periapse distances for each of the orbital motions. The "free" 
terminal parameters are thus argument and time of periapse and longitude of the 
ascending node. The transversal ity conditions associated with the free parameters 
dictate that the initial and final impulses take place at the pericenters of the res- 
pective hyperbolic trajectories, that the pericenters and inclinations of the hyper- 
bolic trajectories be also the pericenters and inclinations of the corresponding 
elliptic orbits. These conditions imply pericenter passage times of t and t^ 
for both ellipse and hyperbola at the earth and target, respectively. 

The trajectory iteration consists of two steps which are alternated until 
convergence is obtained: 

. 

l.(a) Transform R and R to planetocentric coordinates using the 

•X. N 

current values of t^ and t to obtain from the ephemeris the necessary planetary 
velocities : 
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(5.1) 


R 

R 


EX 

TN 


= R 


X 


R 


N 




1. (b) Using the R £X 

straints find R___ and R„ T 
EX TN 


and found in 1. (a) and the terminal con - 

and updated values for t and t such that 

A IM 




and satisfying the conditions that R , R _ v , t define an earth focused Kepler 

JJjX J£X X 

hyperbola with the specified inclination, perigee distance, and perigee passage 
time t Q and that R TN » R TN> t N - define a target focused Kepler hyperbola with 
the specified Inclination, pericenter distance, and pericenter passage time t^. 

The formulas for these calculations are derived below. 

2. (a) Using the updated t . t and the R_„, R_. T found in 1. (b) trans- 

X N EX TN 

form R £x and R^,^ to heliocentric coordinates (again using the ephemeris) 


R x r ex 
r n “ r tn 


+ p 


EX 


+ P, 


TN 


(5.2) 


2. (b) Solve the Lambert problems 

R -* R with transit time t -t 
XI IX 


R T "* R„ with transit time t - t 
IN N I 


to obtain updated values for R and R 


N* 
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3. Return to step 1 and iterate until convergence on R , R , t 

W .X. 

and t are obtained. (12 iterations was typical for 16 digit convergence for 
N 

our test cases). After convergence is obtained, extend the Kepler calculations 

• + 

of 1. (b) to obtain the initial conditions R,_, , R„ for the earth hyperbola and 
' Eo Eo 

final conditions R^, for the target hyperbola. Then the constraints on 

the orbital motions enable us to write 



R 


+ 

Tf 


i 



m e 

< 1+e E> 

2 

Ko 


a E 

U-e E > 


IVf 



(l +e T ) 

1 

2 

R Tf 


a T 

d-e T ) 


l-^Tf 


(5.3) 


The semimajor axes and eccentricities, a^, a^, e^, e^ of the elliptic orbits 
about earth and target are easily calculated from the periapse constraints im- 
posed on these orbits, and the magnitudes of the pericenter velocities (given in 
Eq. (5.3)) are then readily verified. 


The solution of the Lambert problem is routine and we do not give any 

details here. The two point boundary value problem used in step 1. (b) (and also 

(3) 

in step 3) is less well known. Following Pines ' , we wish to fit Kepler hyper- 
bolas to given sphere of influence velocities and pericenter constraints. That is, 
for the earth we are given 


R EX’ l R EX ! r E 


R_ = r , R„ 
Eo o Eo 


inclination = i 


att x 

R * = 0 at t = 0 
Eo 
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while for the target we are given 


letting 


r tn . IR tn I = r T at t= f N 


l R Tff = r f * R xf ' R Tf = 0 att= t f 


V EX R EX R EX’ V TN R TN R TN 


the semimajor axes of these hyperbolas are 


2 

1 2 V EX 



2 

1 _ _2_ V TN 

r T 


(5.4) 


(5.5) 


where we use the subscript H to distinguish hyperbolic elements from elliptic 

elements. Also since r Q and r^, the perigee and target pericenter distances 

are given we can obtain e and e from 

EH TH 


r 

o 


= a 


EH 


(1 - e EH> 


'EH 


a 


EH 


a 


EH 


r 

o 


r f - (i~ e TH ) 


e TH 


a TH ~ r f 
a TH 


(5.6) 


Hence the magnitudes of the angular momentum vectors for these hyperbolas are 
given by 
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h EH = '/“E^H (1 ' e EH> = /Vo^EH* = '.Vo 5 




EH 


h TH = ArV 2 -^* 


TH 


Also, from the definition of angular momentum 


(5.7) 


,2 _ • 2-2 2 ^ . * ,2 

h EH (R EX XR EX ) “ r E V EX " ^EX R EX ) 


h TH (R TN XR TN ) r T V TN " ^TN ' R TN } 


(5.8) 


which yield for R £X • R £X and R TN * R TN 


H • R = / r 2 - h^ 

EX EX v r E EX EH 


R * R = - / r ^ v ^ - h 

TN TN V T TN 


2 

TH 


(5.9) 


when use is made of Eq. (5.7) for h and h„ TI . The difference in the signs 

EH I H 

prefixing the radicals in Eq. (5.8) arises because at t the vehicle exits from 

.2v 

the earth's SOI, while at t it is entering the target's SOI. 

The angular momentum vectors for the hyperbolic orbits must be perpen- 

• • 

dicular to the velocity vectors R EX and R TN « This means that, using the unit 
polar vectors k and k , for earth and target, respectively, we may write 

ili 1 


H EH = °E ^E x R EX> + *E R EX x ^e X R EX> 


H TH = a T x r tn ) + r tn x (k T x r tn ) 


(5.10) 


-22 



where the ct's and /3’s are to be so selected that, using the magnitudes of the 

angular momenta, h and h of Eqs. (5.7) and the specified inclinations 

EH TH 

i and i, 
o f 


“e • H EH “ h EH COS l o *T • H TH = h TH 008 ‘f 


1 1 . * H li 

EH EH EH 


H * H — h 
TH TH TH 


Dropping subscripts temporarily, we thus seek a. and /3 such that 


H = a(kxR) + j3Rx(kxR) 


with 


k • H = h cos i given and H • H = h given. 


From the first of these conditions 


k • H = £(kxR) 2 = /3v 2 sin 2 j = h cos i 


where 


, k • R 

COS 1 

J V 


and hence 


a h cos i 

2.2 

v sin j 


From the second condition 

2 2 2 .2 4 2 2 

a v sin j + /3 v sin j = h 

and, using the expression (5.15) for 0 


(5.11) 

(5.12) 

(5.13) 

(5.14) 

(5.15) 

(5.16) 
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„ , _ . , 2 2 . 
2 1 , 2 h cos i 

a - — — - h - — — 

2 . 2 L 4,4 
v sin j v sm j 


4 . 2. 
v sm j 


-5- fi- 
2 . 2 . 
v sm ] 


, 2 . 
sin j 


This equation for a possesses 


2 2 

no solution for cos i > sin j 


2 2 

one solution, ot = 0, for cos i = sin j 


two solutions, a 


2 . 

cos l 


v|sm j| gin j 


2 2 
for cos i<sin j 


Clearly, the ranges of i for which no solutions exist correspond to inclinations 
which are unattainable from a given velocity asymptote without an additional 
velocity impulse. Although such cases are very real possibilities, we exclude 
them from further consideration here because the special treatment required for 
this inclusion detracts from the primary purpose of the report. 

The angular momentum vector, by definition, is 


and hence 


H =R x R 


RxH = Rv - R (R * R ) 


V V 


But, since H may now be assumed known in terms of a and £ 
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R x H = aR x (k 
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where the dot products may be evaluated from Eqs. (5. 9). 

To update the times t x and t^, we make use of the known semimajor 
axes, eccentricities and perigee passages times (t and t^) of the hyperbolas 
to determine the eccentric anomalies E and E of the SOI points: 

A 1 


r E a EH ( 1_e EH C0Sh E X J r T a TH ( 1-e TH cosh E t ) 


and then use Kepler’s equation to get the times: 


n (t - t ) = e sinh E - E 
E ' X o ; EH XX 


n ft - t ) = e sinh E - E 
T ' N r TH T T 


where 


I 


*E V | ,3 

a EH 


rw 

n„, = / 

T / . _ 


TH 1 


Finally, after the iteration is complete, we carry the Kepler hyperbola calcula- 
tion one step further to obtain 
* + 

R , R by backwards propagation of the converged values for R , 
Eo Eo EX 

**EX’ t X’ t0 t o 


R t j, R by forwards propagation of the converged values for R^, 


R TN’ t0 l f * 


This completes the description of the calculations necessary to implement 

the trajectory iteration, Step I of the double iteration procedure, the updating of 
• • 

t , t , R , R for the next entry into Step I, and the calculation of all para- 

A IN A IN 

meters necessary for Step II of the double iteration procedure. 


(5.26) 


(5.27) 


(5.28) 
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VI. 


PRIMER VECTOR CALCULATIONS 


The ultimate purpose of the optimum three impulse solution is to provide 
a reference trajectory about which to expand a truncated series solution which 
will closely approximate a low thrust trajectory. This approximate trajectory 
should provide an adequate first guess for initiating the iterative numerical so- 
lution of a low thrust trajectory optimization problem using an indirect optimi- 
zation technique. Therefore, even though the three impulse optimization problem 
may be formulated and solved within the context of ordinary calculus, it is help- 
ful to pose the problem in terms of the variational calculus to facilitate imple- 
mentation of the results. 

Past experience with optimization problems involving patched conic tra- 
jectories has indicated that the primer vector and its time derivative are extremely 
sensitive in the proximity of a planet. This presents tremendous difficulties 
in the solution of the boundary value problem because the transversality conditions 
containing these variables evaluated at points near the planets become highly non- 
linear functions of the independent parameters. As a consequence the boundary 

value problem is unstable and virtually Impossible to solve using differential 

( 4 ) 

correction techniques. A method has been found , however, which greatly 
alleviates the sensitivity of the boundary value problem for patched conic tra- 
jectories. This method involves the use of the standard conic equations to re- 
write constraints in the state at a highly sensitive point in terms of the state at 
a less sensitive point. Specifically, it has been found that expressing constraints 
at closest approach to a planet (e.g. , specified passage distance) in terms of con- 
ditions at the sphere of influence is sufficient to permit the solution of the boundary 
value problem. It is for this reason that the optimization problem that follows is 
formulated commencing at exit of Earth's sphere of influence and terminating at 
entrance of the target’s sphere of influence. 
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Proceeding formally as in Reference 4, we define a complete set of 
state variables for the two heliocentric trajectory segments which are joined 
at the midcourse impulse point. In the analysis to follow, the pre -subscripts 
1 and 2 refer to the segments before and after the impulse, respectively. 
Thus, the state equations may be written 



.R 

i 


.R 

i 


i 


V 


( 6 . 1 ) 


for i = 1 and 2. Then, defining r as the time from exit of the Earth's SOI 
to the impulse point, r as the time from the impulse point to entrance of the 

Cl 

target planet's SOI, and s as the independent variable of integration, where 
0 ^ s s i t such that 


2 t= 2 t<0)+ 2 T 8 


1 t (0) + i rs 


one may rewrite the state equations with s as the independent variable of in- 
tegration as follows: 

.V' = - .T .R 

i i .r 3 i 

i 


,R r = .T .V 

l l l 


.t' = .T 
l l 


( 6 . 2 ) 


(6.3) 


r' = 0 
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where the prime denotes differentiation with respect to s. This transformed pro- 
blem fits well within the framework of the indirect method. 


The variational Hamiltonian for the transformed problem is 


h v= I [ t \ • i v ’ + iV i R ' + iV] 


= ) -r\- (Ar * - R ) + A ' -V + A* 

A i L ,r° i V 1 i R 1 it 


Since h is known to be a constant of the motion, it is clear that the bracketed 
v 

term for each leg is also a constant. In fact, the constant 


h= — (A ‘ R) + A • V+ X 
i ,r 3 H V i ; i R l it 


is the more familiar Hamiltonian for the problem in which time is the independent 
variable of integration. The adjoint equations may be written down using the gen- 


eral formula 


A x“- ail v /3X - 


That is, 


A' = - t A 
i V l R 


t A R = t T [ t A V ‘ ( t A V ' i R > l R ] 
i r i r 


i V t ^ 0 


A ’ = - .h 

l T i 
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The equations (6.7) are known to possess analytic solutions. The variables .X 
are obviously constants, the values of which we will subsequently determine to be 
arbitrary. Also, the variables .X ^ are seen to have the solutions 


.X (s) - .X (0) = -.hs. 

IT i T i 

Upon transforming the first two of equations (6. 7) to derivatives with respect 
to time, it becomes clear that 


i 




V 


which leads to the second order equation 


• • 

.A 

i 


V 


UL 


<i A V 


l R) t R 


. JL. 
3 
t r 


, A V 


and since this equation is identical in form to that of the state variational equa- 
tions, its solution will also be of the same form as that of the variational equations. 

In particular, if one partitions the state transition matrix into the four 3x3 matrices 


.A = 


a i R <y 


B = 


w 


C » 


s i R <y 


D = 


s i R <y 


1 Sjiwy ’ 1 a.Rfy ’ 1 ’ f s.R^) 


and if one is given the values of . (t^) and .A v (t^), then the solution to 
(6.10) is 


iVW iVW i A vV 


i A v<v = i c i A v<v + t D i A vV 


This solution is applicable for t less than or greater than t . 


( 6 . 8 ) 


(6.9) 


( 6 . 10 ) 


(6.11) 


( 6 . 12 ) 
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The transversality conditions that are sought are derived from the 
general condition 

2 r 1 

dll+ V [.A ' d.V -.A * d.R + X d.t + .A d .T ] =0 

u V i iV i iti iriJ 

i=l 


in conjunction with the specified constraints on the boundary conditions, where 
II denotes the performance index that is to be minimized. For purposes of 
Illustration, we choose n to be the negative of the ratio of final to initial mass. 
That is 


- A Vvi 


-Av / c 

2 2, 


- Av 3 /0 3. 


n= 1 1 )j[l-(l+k 2 )(l-e 2 2 )] [l-(l+k 3 )(l-e 3 3 )] 


with notation as defined in Section II. For simplicity, we shall assume that the 
launch and capture orbit injection maneuvers are each coplanar with their res- 
pective parking orbits and that the periapses of the hyperbolic and elliptic orbits 
at each terminal are coincident. Under these assumptions the first and third im- 
pulses are written 

Av i * ■Ke+^e'V-v - 

^ v 3 “ “ / ^ X < 1+C T > /,a i( 1_e T ) 


where v aE and v^ are the hyperbolic excess speeds of launch and target hyper- 
bolic trajectories, respectively. The second impulse is equal to the difference in 
the heliocentric velocities on each side of the midcourse impulse, i. e. , 

Av 2 = I 2 V(0)- 1 V(1)| 


(6.13) 


(6.14) 


(6.15) 


(6.16) 


31 - 



Thus, if the launch and target parking orbits are specified in terms of semi-major 
axis and eccentricity, it is seen that II may be written functionally as 


n= n ( v w g> /(i). 


2 V <°>- W 


(6.17) 


such that 


dn= " V «E dV “E + \ T + V (1 > ' dlV<1) + V«»‘ V<°> <6 - 18) 


with 


n 


an 

v „ BAv 


00 E 


» E V 1 + 2M E /a E (l-e E ) 


n 


Bn 


00 T 


v m B Av 
«x 


n 2 V(0) n l v(D 


3 /vJ T +2^ T /a T (l-e T ) 

2 V (°)-iV(l) , n 


Av„ BAv„ 


(6.19) 


where 


-Av./ c. 


*17 r -av/'' -Av /c (1+k.) -Av./c 

5^- - [l-a-Ht >a-e 


(6.20) 


i, j, m = 1, 2, 3; i / j / m 


The specified boundary conditions of the problem affect the transversality 
conditions through the differentials of the state appearing in equation (6.13). Em- 
ploying the notation of Section II, one will recognize that 
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( 6 . 21 ) 


i R(0) r x p EX + p EX ; 1 V <°) R x r ex + p ex 


2 R(!) R n R tn + P tn ’ 2 V(1) R N R TN + P TN 


and since the planetary positions and velocities may be considered functions only 
of time, we may write 


d jBfO) - dR EX + P Ex d lt (0) ; d lV (0) = dR EX + P^MP) 


d 2 R(1) “ dR TN + P TN d 2 t(1) ! d 2 V(1) =dR TN + P T N d 2 t(1) 


( 6 . 22 ) 


A standard formula for decomposing the differential of a position or velocity 
vector into the differentials of polar coordinates is given in Reference 4, and 
leads to 


, kg x h £ 

dR EX =(£ E XR EX )dn E + (fi E xR BX )d "EX + (“inlT XR Ex) dl 

E 


E 


R 


dR„ = 


EX 


k xh 


EX v 


EX 


d V EX +(k E X R EX } d °E + ^E X R EX> (< * W EX~ dy EX )+ ( sin l X R Ex) d l E 

E 


k T x h T 


dR TN X R TN )d ^T + ^T X R TN )d W TN + ( sin i X R TN ) d 1 


(6.23) 


R 


TN 


dR TN = d V TN + x R TN> d °r + ‘“t X R TN ; <d “TN ' d ’'TN* 


/ k T x • \ 

+ l-itoT7 XR TN) di T 

T 


- 33 - 



where O, 60, i, and y are the longitude of ascending node on the equator, the 
argument of position at the sphere of influence, the inclination relative to the equa- 
tor, and the flight path angle at the sphere of influence, respectively, of the plane- 
tocentric hyperbola, and h denotes a unit vector in the direction of the planeto- 
centric angular momentum. Boundary conditions which link the two trajectory 
segments are 

2 R(0) « x Ra); 2 t(0) = 1 t(i) 


which imply 


d 2 R(0) = d 1 R(1) ; d 2 t(0)=d 1 t(l). 

Because the heliocentric flight times r are parameters 

d.r(l) = d ,r(0) = d t 

and these differential may be replaced in (6. 13) using the relations 

d.r = d.t(l) -d.t(0) 

The times at which the spheres of influence are crossed are generally of little 
immediate importance to the trajectory analyst. Of greater interest are the times 
at which the launch and target injection maneuvers are performed, and it is useful 
to replace the differentials of t(0) and t(l) with the differentials of the launch 

X u 

time t and the target injection time t,. Denoting t and t as the times 

O f 35 E 00 x 

within the spheres of influence of Earth and the target planet, respectively, then 

l^-V^E > 2 tW “ t f- t -T 
d 1 t(0)=dt o+ dt wE ; d 2 t(l)=dt f -d te#T 


(6.24) 


(6. 25) 


(6.26) 


(6.27) 


(6.28) 
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where 


*.E = (U E /V .E >< e EH Slnh W 


e EH " 1 + a E (1 ' e E ) V «E /W E 


f E '°"h‘ 1 [(l^ E ''. 2 E /(i E )/e EH ] 


t ®T ^ U T^ V ®T^ C TH K ' n ^ ~ 


e TH = 1+a T (1 ' e T> V -T /U T 


f = cosh * fil+r v 2 /u )/e,„,, 1 
T L' T «T “ T THJ 


Thus the times within the spheres of influence are seen to be functions of the hyper- 
bolic excess speeds plus other known functions and we may write 

at _ 

dt = . ■ d v 

«°E ov „ ®E 
«E 


St T 

dt -T “ dv -T 


where 


dt «E 3t «E 2 

dv„ 'v/ 2 . . . L'~EH ““ *E * E '~EH E' 

®E 00 E v #E e EH smhf E 


[< e EH°° eh f E _1) ; E- (e EH- C ° 8hf E )(1 - e E )a E ] 


dt »T _ 3t «T ^_ 2 

dV «T V “T v„ 2 t e TH stnh t T 


[ (e TH C ° Shf T* 1) ^T " (e T H " COShf T )(1 " e T) a T ] 
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(6.29) 


(6.30) 


(6.31) 



Finally, it is not difficult to note that the speed and the flight path angle at the 
sphere of influence are functions only of the hyperbolic excess speed and other 
known parameters. Consequently, one may replace the differentials of these two 
functions, which appear in (6.23), with the differential of the excess speed. 
Employing the relations 


r 

o 



< 1 ~ e E ); 

(l-e T ); 


2 2 

V = v 
oH 00 E 

2 2 

V = V 

fH »T 


+ 2u /r 
^E o 

+ 2/Lt T /r f 


(6.32) 


r E V EX ° 0S y EX = VoH ; 'T V TN C ° S y TN = r f V fH 


one obtains 


2 2 

'EX , A . V »E L V oH _V EX„ 

d v_ - (h„ x R pv )dy T?v = ~ LR~~ + • R 


v £X EX ' E 


EX' ~~ ' EX 2 
V oH 


EX 


] dv . 


(R , R ) EX -I »E 
n EX EX ; 


(6.33) 


R 


TN 


TN 


dV TN " X R TN* dy TN 


00 X 


fH 


t R ™ + (R, 


2 2 
V -V 

fH TN 


R 


TN R TN ) 


TN 


] 


d v 


Upon substituting the results of equations (6.18) - (6.33) into the general 
equation (6.13) and setting the terms associated with the remaining independent 
differentials to zero individually yields the appropriate transversality conditions 
for our problem. 
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(a) 

<b) 

(c) 

(d) 

(e) 

(f) 

(g) 

(h) 

(D 


( 2 A V ( 0 ) - 1 A V (1) )' d 2 R(0) '° 

( n ,ve»'2 A v <0 ) )- d 2 v<0) = 0 

( 2 A V (0) * 2^ (0) ' 1 A V (1) * !V(1)- 2 A v (0) * 2 V(0) + i A v ( 1) • 1 V(l))d 2 t(0) = 0 

(l A V (0) -“EX-lV 0 * • R Ex) dt o + ( 2 A V (1 » • R TN-2 A V (1) - R T N ) dt f = 0 
V ( R EX X lV 0, ' R EX X l A V <0) ) dn E = 0 
fi E ’ ( R EX X 1 A V (0) ' R EX X i A V<°») d “EX = ° 


k x h 
E E 

sin 


( R EX X lV 0 ) - R EX X lV 0) ) di E = ° 


L v 


2 2 

V v - V 

•E r« . , . OH EX 


00 E v 


_ R EX * lV 0) + 


oH 


^EX * R EX ) 


( r ex • lV»)] 


(i) 

(k) 

a) 


dt 

dv 


00 E 
®E 


(l A V<°>' R EX-l A V ( 0 ) - R Ex)} d ''.E = ° 


*T • ( R TN X 2 V - R TN X 2 A V< 1 >) * °T = ° 


V ( R TN X 2 A V (1 »- R TN X 2V 1 >) 


d co = 0 
TN 


k T X h T 
sin L 


( R TN X 2 A V (1) - R TN X 2V 1 >) di T = ° 
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(6.34) 



<m) {n 


+ °°T 

v _ 2 

ao T v 

fH 


[*TN ' 2V I + 


2 2 
V fH ~ V TN 

r tn* r tn 


( R Tn' 2 A V (1) )] 


dt 

d 


v (2^V (1) ’ R TN " 2 A V (1) ’ R TN)) dV ®T 
coT 


= 0 


(6.34) 

cont 


Applying the standard interpretation to these equations, we regard a boundary con- 
dition to imply that the differential of the state variable which is specified must be 
zero, thereby satisfying identically the associated equation in (6.34). If, however, 
no constraint is imposed, then the differential is arbitrary and its coefficient must 
vanish, and it is this vanishing of the coefficient that represents the transversality 
condition. For the problem under consideration, we assume that no constraints are 
placed on the time and location of the midcourse impulse nor on the velocities on 
each side of the impulse. Consequently, the coefficients of the differentials in Equa- 
tions (a) - (d) in (6.34) must each vanish individually. From (a) we determine 
that the time derivative of the primer is continuous over the second impulse. The 
equations (b) and (c) define the primer vector explicitly at the midcourse impulse. 
They indicate that, for the performance Index of our problem, the primer is con- 
tinuous over the impulse and directed along the impulse. The use of the results of 
(a) - (c) in equation (d) leads to the conclusion that, on the optimal trajectory, 
the derivative of the primer at the second impulse is perpendicular to the impulse 
and therefore to the primer itself. The equations (e) - (m) yield transversality 
conditions associated with parameters at the launch and target planets, and are 
used as described below. Note that X has disappeared entirely from the equations; 
hence the earlier statement that X is arbitrary. 

The evaluation of the adjoint variables along the trajectory is accomplished 
with the use of equations (6. 12) which require that both the primer vector and its 
time derivative be known at one instant in time. Parts (b) and (c) of (6.34) 
give the primer vector at one point in time for each trajectory segment. To define 
the time derivative of the primer vector on each segment at the same time requires 
a total of six independent equations (i. e. , three for each segment). For this purpose 
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we choose parts (f), (g) and (i) for use on the first trajectory segment and 
parts (j), (k) and (m) for the second segment. The equations (h) and (1) 
are not used for this purpose because the trajectory iteration procedure des- 
cribed in Section V assumes that inclination at both ends are specified. In de- 
fining the partitioned transition matrices for the first segment, let the t of 
equations (6.11) represent the time at the impulse and let t represent the 
time at exit of the Earth’s sphere of influence. Then using (6.12) to substitute 

o 

for A v (0) and v (0) in (f), (g) and (i), and us ing the vector operation 

identity 

X ’ (YxZ) = (X x Y) * Z 


yields the equations 

<*E * R EX> ' <l°l A v a > + * R EX> ' <1 A 1 A V (1) + A A V (1 » “ ° 

x r ex ) ' <i c i A v (1) + i D i A v (1)) - ^E x R EX> ' < 1 A 1 A V (1) + ,V,«! “ 0 


V 

~ r ®e /• 

n v L~2 ( R EX + 

" E V oH 


2 2 
v - v 
oH EX 

R * R 
EX EX 


r ex)-7T7 r ex]' (iVv^lV 
®E 


dt 


® E 


dv 


00 E 


r ex- (i c iV 1) + i d i A v (1) ) ■ 0 


Upon noting the identity 


Y • (QX) = (Q Y) * X 


where Q is a matrix and X and Y are vectors, and reorganizing equations 
(6.35), one obtains 


(6.35) 

v (1 >) 


(6.36) 
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DT ^E X V , -l BT ^ xk EX ) ]- 1 A V< 1 > = [l AT < i; E Xii EX , -l CT < i; E xR EX ) }A< 1) 


r _T 

Li 


T- 


.T - 


Ll D (h E XR EX ) 'l B (h E XR EX ) J- 1 A v a) = [ 1 A T (h E xi ?Ex )- i C T (h E xR EX )]- 1 A v ^ 


f -E d t r + b T r °° E (■ 

Idv _ 1 D R EX 1 B L 2 V 
“ E V oH 


r ex + 


2 2 
V oH ~ V EX 

R EX r EX 


dt 


) 00 TT • • ”*p ■ * 

" dV~ R EX-if 'l A V (1) = n v 


“ E 


c- R „....A T r^(R EX , T °» . :> r 


dt 

" Idv „ 1 C “EX ' 1“ L 2 


00 E 
dt 


(6.37) 
00 E 


oH 


p . p EX/ dv _ 

r ex r ex 00 e 


\ “E .. -H 

)‘ 5 r: R Ex]}'iV l) 


But these three equations, when written in matrix form, are 


Q 1 A y (l) = Y 


(6.38) 


which has the solution 


1 A y (l) = Q -1 Y 


(6.39) 


Proceeding similarly for the second leg with t^ again representing the time at 
the impulse but t representing the time at entry of the target sphere of influence, 
the equations for the second segment corresponding to (6.37) are 


[ 2 D T (k T xR TN )- 2 B T (k T xR TN )]- 2 A v (0) = [ 2 A T (k T x RTN )- 2 C T (k T x RTN )]. 2 A v (0) 


[ 2 D X XR TN> A T <4 XR TnO' 2 A V (0> l 2 AT < ii T XR TN» -2 CT < i T XR TN ) ]- 2 A V (0) 


{ 


-°^00T T • T r V coT 

dv 2 TN 2 L 2 


00 X 


V — V 

fH TN 


fit + 

^ TN R * R 


dt 

^ . “T 




R TN r^r- R TN J/° 2 A v (0)=n ’ 


fH 


TN TN 


coT 


(6.40) 


oo T 


dt 


r °°t 
tdv m 

00 T 


oo x T 


t r »t 


2 TN 2 


C'R^-jr j — — (A+ - R, 

L 2 \ TN n • r “ 

v, T\r tm 00 T 

fH 


2 2 
V -V 

fH TN 


. dt m 
\ ooT - 
)+- H 


Rn 


TN R • r TN J dv _ TN 
TN TN “ T 


JKV* 
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The procedure for evaluating the primer vector and its time derivative 
along a trajectory resulting from the iteration of Section V is now clear. From 
equations (6.34) parts (b) and (c), the primer vector at the impulse point is 
evaluated. This, in conjunction with information pertaining to the trajectory 
(including the transition matrices for the two segments), is used in equations 
(6.37) and (6.40) to obtain the A^'s on both sides of the impulse. These will, 
in general, not satisfy either (6.34), part (a) or part (d). A differential 
correction scheme is then employed to adjust the time and position of the second 
impulse until the four equations represented by parts (a) and (d) are satisfied. 

The preceding is sufficient to cover the problem with specified launch 
and target arrival dates and specified inclination at the two terminals. One may, 
however, optimize any combination of these parameters using the remaining equa- 
tions in (6.34) by simply expanding the order of the boundary value problem im- 
plied in the preceding paragraph. For example, the planetocentric inclinations at 
launch and/or target may be optimized by driving to zero the coefficients of (h) 
and (1), respectively. One simply adds these coefficients to the other list of 

boundary conditions, parts (a) and (d), and treats i and i as independent 

£4 I 

parameters along with the position and time of the impulse. Similarly, the co- 
efficients of dt and d t . in part (e) must vanish independently if t and t, 
o t 01 

are to be optimized. If both t Q and t^ are left open but the flight time t^ - t Q 

is constrained, then dt = dt„, and the sum of the two coefficients in (e) must 

o i 

vanish, thereby providing a single boundary condition with either t Q or t^ serving 
as the independent parameter. 

Once the boundary value problem is solved, one may elect to obtain the 
primer vector within the planetocentric segments of the trajectory. This may be 
accomplished by proceeding in both directions from the second impulse along the 
heliocentric arcs using equations (6. 12). At the spheres of influence the dis- 
continuities in the Ay's are applied as defined by equations (2.7). Then the 
equations (6.12) once again apply for propagating the adjoint variables along the 
planetocentric conic trajectory. 
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VII. NUMERICAL RESULTS 


To test the approach described in the preceding sections, an Earth-Mars 
transfer was chosen for illustrative purposes. A trajectory leaving Earth orbit 
at noon, January 21, 1977 and arriving in Mars orbit 270 days later on October 
18, 1977 was selected. For computational purposes, the altitudes of the initial 
and final circular orbits were both taken to be zero. The assumed characteristics 
of the first stage included a jet exhaust speed c^ of 4.4 km/sec and a structural 
factor k^ of 0.1. The corresponding parameters of the second and third stages 
were assumed to be as follows : 

c = 30 km/sec, k =0.03, c = 3 km/sec, k =0.1. 

Z Z o o 

The heliocentric ecliptic transfer angle for this mission is approximately 
315 degrees, commencing at a longitude of 121 degrees and reaching Mars at a 
longitude of 76 degrees. The longitude of ascending node of Mars orbit on the 
ecliptic plane is about 49 degrees while the longitude of aphelion of the heliocen- 
tric ecliptic two-body trajectory connecting Earth and Mars occurs at a longitude 
of 65 degrees. Thus, it is clear that both apsides and both nodes are contained 
within the heliocentric ecliptic transfer and therefore each of the four programmed 
starting guesses for tj and R^ may be invoked for this case. But, because the 
longitudes of the line of nodes and the line of apsides differ by only 17 degrees, 
the starting guesses obtained by using the ascending node and the apoapse point 
are very close as are the guesses associated with the descending node and the 
periapse point. 

The solution to the problem posed may be characterized as follows. 

About 99 days following departure from Earth the intermediate impulse of 3.8 
km/sec is applied at a longitude of 259 degrees and a latitude of -1 degree. 

The heliocentric distance at the impulse is about 0.64 AU. The impulses at 
launch and arrival are 4.36 km/sec and 4.16 km/sec, respectively, and the 
final to initial mass ratio is 0.047. 
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The programmed starting guess closest to this solution is the periapse 
point which lies 14 degrees in longitude from the solution. The periapse dis- 
tance of the heliocentric ecliptic trajectory is about 0.425 AU. Although this 
would seem to be a reasonably close first guess, it was found that the program 
would not converge to the solution from the ecliptic periapse point, nor from 
any of the other programmed first guesses for that matter. A principal reason 
for this is that the mathematical definition of the performance index, as stated 
in equation (6.14), is not always descriptive of the actual situation. For example, 
it is clear that if a velocity impulse is sufficiently large, the propellant mass 
requirements can exceed the difference in the initial stage mass and the engine 
and tankage mass such that the corresponding term within the square brackets of 
(6.14) becomes negative. Although this represents a physically unrealizable 
situation, the circumstance should not hamper convergence because the partial 
derivatives of the performance index with respect to the velocity impulse will 
still be negative. However, if two velocity impulses are sufficiently large that 
the corresponding bracketed terms in (6.14) are both negative, then the partial 
derivatives of II with respect to each of those impulses will be positive and 
hence indicate that the desired solution is in the direction of increasing velocity 
impulses. The result is that subsequent iterations are unable to move the trial 
impulse point and time to a region where physically meaningful solutions are 
available, and the technique simply fails to converge. The importance of this 
problem is placed in perspective when one realizes that each of the four pro- 
grammed starting guesses for the problem posed results in precisely this situa- 
tion. 

This problem may be circumvented by first solving the problem with the 
structural factors, k , set to zero. This prevents the individual stage payload 
fractions from becoming negative and thereby avoids the problem. The result- 
ing solution should then represent a reasonable first guess to the problem with 
non-zero structural factors. This approach was attempted with each of the 


- 43 - 



programmed starting guesses. As might be expected for this case, starting 
from either the ascending node or the aphelion point, the program failed to con- 
verge because the starting guesses were far from the known solution. Starting 
from either the descending node or the periapse point, successive iterations 
did proceed toward the solution. However, as the transfer angle on the second 
trajectory segment (i.e., following the intermediate impulse) approached 180 
degrees (note that either first guess results in an initial angle of greater than 
180 degrees) the improvement per iteration slowed significantly and finally 
halted altogether. This was due to the fact that the near 180 degree three-dimen- 
sional transfer between specified points represents a highly sensitive region in 
which the behavior of the state and adjoint variables are transient. This results 
in irregularities in the multi-dimensional surface representing the end conditions, 
and these irregularities can present an insurmountable obstacle using techniques 
which depend on local slopes of that surface, such as the one presented here. 

Various approaches were attempted to alleviate this latter problem. Direct 
minimization of the performance index (a feature of the iterator employed ^) 
starting from the periapse point was unsuccessful because progress halted at the 
180 degree transfer point. Employing a first guess on the heliocentric ecliptic 
trajectory that was in the appropriate region in terms of angular position was also 
unsuccessful because the curvature of the surface is such that the iterator over 
corrects in longitude causing the next trial trajectory to cross the offensive bar- 
rier once again. 

It is believed that the specific problem investigated displays the major 
difficulties that one would expect to encounter with the approach described herein. 
The results of this case imply that, to use the approach effectively, a better method 
for obtaining initial guesses for the time and location of the midcourse impulse 
must be developed, because the radius of convergence for some problems will 
obviously be very small. Alternatively, improved techniques for controlling the 
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magnitude and direction of the steps taken by the iterator could be developed. It 
is suggested that more extensive use of the program be made for a variety of 
missions as the need and opportunity arises, and that results of these studies be 
employed to assist in the development of the necessary improved starting techniques. 
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